In [1]:
! pip install geopandas
! pip install matplotlib
! pip install matplotlib-scalebar
! pip install osmnx
Requirement already satisfied: geopandas in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (0.14.3)
Requirement already satisfied: fiona>=1.8.21 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from geopandas) (1.9.6)
Requirement already satisfied: packaging in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from geopandas) (23.1)
Requirement already satisfied: pandas>=1.4.0 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from geopandas) (2.2.2)
Requirement already satisfied: pyproj>=3.3.0 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from geopandas) (3.6.1)
Requirement already satisfied: shapely>=1.8.0 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from geopandas) (2.0.4)
Requirement already satisfied: attrs>=19.2.0 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas) (23.1.0)
Requirement already satisfied: certifi in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas) (2023.7.22)
Requirement already satisfied: click~=8.0 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas) (8.1.7)
Requirement already satisfied: click-plugins>=1.0 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas) (1.1.1)
Requirement already satisfied: cligj>=0.5 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas) (0.7.2)
Requirement already satisfied: six in /Users/ezrabanks/Library/Python/3.11/lib/python/site-packages (from fiona>=1.8.21->geopandas) (1.16.0)
Requirement already satisfied: numpy>=1.23.2 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from pandas>=1.4.0->geopandas) (1.25.2)
Requirement already satisfied: python-dateutil>=2.8.2 in /Users/ezrabanks/Library/Python/3.11/lib/python/site-packages (from pandas>=1.4.0->geopandas) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from pandas>=1.4.0->geopandas) (2023.3)
Requirement already satisfied: tzdata>=2022.7 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from pandas>=1.4.0->geopandas) (2024.1)
Requirement already satisfied: matplotlib in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (3.8.4)
Requirement already satisfied: contourpy>=1.0.1 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib) (1.2.1)
Requirement already satisfied: cycler>=0.10 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib) (4.51.0)
Requirement already satisfied: kiwisolver>=1.3.1 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib) (1.4.5)
Requirement already satisfied: numpy>=1.21 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib) (1.25.2)
Requirement already satisfied: packaging>=20.0 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib) (23.1)
Requirement already satisfied: pillow>=8 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib) (10.3.0)
Requirement already satisfied: pyparsing>=2.3.1 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib) (3.1.2)
Requirement already satisfied: python-dateutil>=2.7 in /Users/ezrabanks/Library/Python/3.11/lib/python/site-packages (from matplotlib) (2.8.2)
Requirement already satisfied: six>=1.5 in /Users/ezrabanks/Library/Python/3.11/lib/python/site-packages (from python-dateutil>=2.7->matplotlib) (1.16.0)
Requirement already satisfied: matplotlib-scalebar in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (0.8.1)
Requirement already satisfied: matplotlib in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib-scalebar) (3.8.4)
Requirement already satisfied: contourpy>=1.0.1 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (1.2.1)
Requirement already satisfied: cycler>=0.10 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (4.51.0)
Requirement already satisfied: kiwisolver>=1.3.1 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (1.4.5)
Requirement already satisfied: numpy>=1.21 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (1.25.2)
Requirement already satisfied: packaging>=20.0 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (23.1)
Requirement already satisfied: pillow>=8 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (10.3.0)
Requirement already satisfied: pyparsing>=2.3.1 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (3.1.2)
Requirement already satisfied: python-dateutil>=2.7 in /Users/ezrabanks/Library/Python/3.11/lib/python/site-packages (from matplotlib->matplotlib-scalebar) (2.8.2)
Requirement already satisfied: six>=1.5 in /Users/ezrabanks/Library/Python/3.11/lib/python/site-packages (from python-dateutil>=2.7->matplotlib->matplotlib-scalebar) (1.16.0)
Requirement already satisfied: osmnx in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (1.9.2)
Requirement already satisfied: geopandas>=0.12 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from osmnx) (0.14.3)
Requirement already satisfied: networkx>=2.5 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from osmnx) (3.3)
Requirement already satisfied: numpy>=1.20 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from osmnx) (1.25.2)
Requirement already satisfied: pandas>=1.1 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from osmnx) (2.2.2)
Requirement already satisfied: requests>=2.27 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from osmnx) (2.31.0)
Requirement already satisfied: shapely>=2.0 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from osmnx) (2.0.4)
Requirement already satisfied: fiona>=1.8.21 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from geopandas>=0.12->osmnx) (1.9.6)
Requirement already satisfied: packaging in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from geopandas>=0.12->osmnx) (23.1)
Requirement already satisfied: pyproj>=3.3.0 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from geopandas>=0.12->osmnx) (3.6.1)
Requirement already satisfied: python-dateutil>=2.8.2 in /Users/ezrabanks/Library/Python/3.11/lib/python/site-packages (from pandas>=1.1->osmnx) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from pandas>=1.1->osmnx) (2023.3)
Requirement already satisfied: tzdata>=2022.7 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from pandas>=1.1->osmnx) (2024.1)
Requirement already satisfied: charset-normalizer<4,>=2 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from requests>=2.27->osmnx) (3.2.0)
Requirement already satisfied: idna<4,>=2.5 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from requests>=2.27->osmnx) (3.4)
Requirement already satisfied: urllib3<3,>=1.21.1 in /Users/ezrabanks/Library/Python/3.11/lib/python/site-packages (from requests>=2.27->osmnx) (1.26.16)
Requirement already satisfied: certifi>=2017.4.17 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from requests>=2.27->osmnx) (2023.7.22)
Requirement already satisfied: attrs>=19.2.0 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas>=0.12->osmnx) (23.1.0)
Requirement already satisfied: click~=8.0 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas>=0.12->osmnx) (8.1.7)
Requirement already satisfied: click-plugins>=1.0 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas>=0.12->osmnx) (1.1.1)
Requirement already satisfied: cligj>=0.5 in /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas>=0.12->osmnx) (0.7.2)
Requirement already satisfied: six in /Users/ezrabanks/Library/Python/3.11/lib/python/site-packages (from fiona>=1.8.21->geopandas>=0.12->osmnx) (1.16.0)
In [13]:
pip install matplotlib-scalebar
Collecting matplotlib-scalebar
  Obtaining dependency information for matplotlib-scalebar from https://files.pythonhosted.org/packages/a9/9e/22930e3deb2c374f47c6633aff9f6f379f8c421ab868fff3b4f85eac8b8a/matplotlib_scalebar-0.8.1-py2.py3-none-any.whl.metadata
  Downloading matplotlib_scalebar-0.8.1-py2.py3-none-any.whl.metadata (13 kB)
Requirement already satisfied: matplotlib in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from matplotlib-scalebar) (3.7.1)
Requirement already satisfied: contourpy>=1.0.1 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (1.0.5)
Requirement already satisfied: cycler>=0.10 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (0.11.0)
Requirement already satisfied: fonttools>=4.22.0 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (4.25.0)
Requirement already satisfied: kiwisolver>=1.0.1 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (1.4.4)
Requirement already satisfied: numpy>=1.20 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (1.24.3)
Requirement already satisfied: packaging>=20.0 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (23.0)
Requirement already satisfied: pillow>=6.2.0 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (10.2.0)
Requirement already satisfied: pyparsing>=2.3.1 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (3.0.9)
Requirement already satisfied: python-dateutil>=2.7 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from matplotlib->matplotlib-scalebar) (2.8.2)
Requirement already satisfied: six>=1.5 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from python-dateutil>=2.7->matplotlib->matplotlib-scalebar) (1.16.0)
Downloading matplotlib_scalebar-0.8.1-py2.py3-none-any.whl (17 kB)
Installing collected packages: matplotlib-scalebar
Successfully installed matplotlib-scalebar-0.8.1
Note: you may need to restart the kernel to use updated packages.
In [15]:
pip install osmnx
Collecting osmnx
  Obtaining dependency information for osmnx from https://files.pythonhosted.org/packages/55/a7/cf80815256c78878dc0749496652ba5e1d54445b6a01951f6a6f345bcc3f/osmnx-1.9.2-py3-none-any.whl.metadata
  Downloading osmnx-1.9.2-py3-none-any.whl.metadata (4.9 kB)
Requirement already satisfied: geopandas>=0.12 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from osmnx) (0.14.2)
Requirement already satisfied: networkx>=2.5 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from osmnx) (3.1)
Requirement already satisfied: numpy>=1.20 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from osmnx) (1.24.3)
Requirement already satisfied: pandas>=1.1 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from osmnx) (1.5.3)
Requirement already satisfied: requests>=2.27 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from osmnx) (2.31.0)
Requirement already satisfied: shapely>=2.0 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from osmnx) (2.0.1)
Requirement already satisfied: fiona>=1.8.21 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from geopandas>=0.12->osmnx) (1.9.5)
Requirement already satisfied: packaging in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from geopandas>=0.12->osmnx) (23.0)
Requirement already satisfied: pyproj>=3.3.0 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from geopandas>=0.12->osmnx) (3.6.1)
Requirement already satisfied: python-dateutil>=2.8.1 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from pandas>=1.1->osmnx) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from pandas>=1.1->osmnx) (2022.7)
Requirement already satisfied: charset-normalizer<4,>=2 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from requests>=2.27->osmnx) (2.0.4)
Requirement already satisfied: idna<4,>=2.5 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from requests>=2.27->osmnx) (3.4)
Requirement already satisfied: urllib3<3,>=1.21.1 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from requests>=2.27->osmnx) (1.26.16)
Requirement already satisfied: certifi>=2017.4.17 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from requests>=2.27->osmnx) (2024.2.2)
Requirement already satisfied: attrs>=19.2.0 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas>=0.12->osmnx) (22.1.0)
Requirement already satisfied: click~=8.0 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas>=0.12->osmnx) (8.0.4)
Requirement already satisfied: click-plugins>=1.0 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas>=0.12->osmnx) (1.1.1)
Requirement already satisfied: cligj>=0.5 in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas>=0.12->osmnx) (0.7.2)
Requirement already satisfied: six in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas>=0.12->osmnx) (1.16.0)
Requirement already satisfied: setuptools in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas>=0.12->osmnx) (68.0.0)
Downloading osmnx-1.9.2-py3-none-any.whl (107 kB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 107.4/107.4 kB 2.5 MB/s eta 0:00:00a 0:00:01
Installing collected packages: osmnx
Successfully installed osmnx-1.9.2
Note: you may need to restart the kernel to use updated packages.
In [17]:
import geopandas as gpd
import pandas as pd
import numpy as np

import matplotlib.pyplot as plt
from matplotlib_scalebar.scalebar import ScaleBar

import osmnx as ox

# Set display option to show all columns
pd.set_option('display.max_columns', None)

Lab 1: Geographic Data and Spatial Models¶

Welcome to today's lab session!

In this session, we will delve into the functionalities of two powerful Python libraries: Geopandas and GurobiPy.

  1. Geopandas: Geopandas is a Python library that facilitates working with geospatial data. It extends the capabilities of Pandas, a widely used data manipulation library, to support geographic data types and operations. Throughout this lab, we will learn how to read, write, visualize, and analyze geographic data using Geopandas.

  2. GurobiPy: GurobiPy, or Gurobi Optimization Python Interface, is a Python package that provides access to the Gurobi Optimizer, a powerful mathematical optimization solver. Gurobi is widely used in various industries for solving linear programming, integer programming, quadratic programming, and other optimization problems. In this lab, we will explore the basics of using GurobiPy to formulate and solve optimization problems.

  3. Knapsack Problem: As part of our exploration, we will delve into a classic optimization problem known as the Knapsack Problem. This problem involves selecting a combination of items to maximize value while staying within a given weight constraint. We will apply the concepts learned from Geopandas and GurobiPy to solve practical scenarios, such as selecting parcels to buy under various conditions, resembling real-world decision-making scenarios.

Throughout the lab, we encourage you to actively engage in hands-on exercises, experiment with code snippets, and explore additional functionalities of Geopandas and GurobiPy. Feel free to ask questions, collaborate with your peers, and enjoy the journey of discovering the capabilities of these powerful Python libraries.

Lab Structure¶

Lab 1 is due in two weeks and consists of two chapters.
Chapter 1 is to figure out which parcel you'd like to buy within given budget.
Chpater 2 is to figure out the distribution of the distance between wildfire ignitions and major roads.

Chapter 1. Which parcel do you want to buy?¶

Chapter 1 performs an exploratory data analysis on Parcel Data through Geopandas.
Then, utilizing given information, a Knapsack problem will be solved using GurobiPy solver.
You'll get to know which parcels you'd like to buy with given budget.

Step 1. Explore Santa Barbara Parcel Data¶

In [8]:
# Read parcel data in
parcels = gpd.read_file("Lab1Data/parcel_SB.shp")

# Plot your parcel data
FIGSIZE = (15,5)
In [9]:
# Check the data types of your parcel data
parcels
Out[9]:
OBJECTID APN LAYER Situs1 Situs2 Acreage LandUse UseCode PctTransf LandValue StrImpr LivImpr NetSecVal Net_Impr YearBuilt Bedrooms Bathrooms geometry
0 131101 041-350-024 Ground 2315 EDGEWATER WAY SANTA BARBARA, CA 93109 0.52 SINGLE FAMILY RESIDENCE 0100 100.0 288614 331840 0 620454 331840 1967 4 3.00 POLYGON ((6039997.884 1971525.657, 6039982.106...
1 131102 041-350-012 Ground 2307 EDGEWATER WAY SANTA BARBARA, CA 93109 0.18 SINGLE FAMILY RESIDENCE 0100 26.6 550902 1085141 0 1636043 1085141 2006 3 3.50 POLYGON ((6040044.305 1971477.086, 6040032.277...
2 131103 041-350-015 Ground 2211 EDGEWATER WAY SANTA BARBARA, CA 93109 0.59 SINGLE FAMILY RESIDENCE 0100 100.0 1641000 175000 0 1816000 175000 1954 2 2.50 POLYGON ((6040237.820 1971288.705, 6040179.849...
3 131104 041-350-016 Ground 2201 EDGEWATER WAY SANTA BARBARA, CA 93109 0.50 SINGLE FAMILY RESIDENCE 0100 100.0 1644611 1482404 0 3127015 1482404 1948 3 3.00 POLYGON ((6040320.008 1971271.353, 6040268.735...
4 131105 041-350-017 Ground 2111 EDGEWATER WAY SANTA BARBARA, CA 93109 0.76 SINGLE FAMILY RESIDENCE 0100 0.0 63587 156682 0 213269 156682 1950 4 3.00 POLYGON ((6040450.139 1971243.880, 6040398.978...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
1867 250105 045-290-008 Condo Third Floor 30 BARRANCA AVE 3 SANTA BARBARA, CA 93109 0.00 CONDOS,COMMUNITY APT PROJS 0323 100.0 46816 70592 0 117408 70592 1973 2 1.75 POLYGON ((6048328.471 1972754.300, 6048325.112...
1868 250106 045-290-009 Condo Third Floor 36 BARRANCA AVE 6 SANTA BARBARA, CA 93109 0.00 CONDOS,COMMUNITY APT PROJS 0323 100.0 750000 668000 0 1418000 668000 1973 2 1.75 POLYGON ((6048287.470 1972741.338, 6048291.580...
1869 250107 045-290-010 Condo Third Floor 36 BARRANCA AVE 5 SANTA BARBARA, CA 93109 0.00 CONDOS,COMMUNITY APT PROJS 0323 100.0 309553 206366 0 515919 206366 1973 2 2.00 POLYGON ((6048283.551 1972700.429, 6048271.984...
1870 250108 045-290-011 Condo Third Floor 40 BARRANCA AVE 3 SANTA BARBARA, CA 93109 0.00 CONDOS,COMMUNITY APT PROJS 0323 100.0 124884 187149 0 312033 187149 1973 2 1.75 POLYGON ((6048234.146 1972677.397, 6048238.386...
1871 250109 045-330-008 Condo Third Floor 101 OCEANO AVE 5 SANTA BARBARA, CA 93109 0.00 CONDOS,COMMUNITY APT PROJS 0332 100.0 305105 298615 0 603720 298615 1961 2 1.75 POLYGON ((6048280.057 1972978.876, 6048296.690...

1872 rows × 18 columns

In [10]:
# Check the Coordiante System of your parcel data
parcels.crs
Out[10]:
<Projected CRS: EPSG:2229>
Name: NAD83 / California zone 5 (ftUS)
Axis Info [cartesian]:
- X[east]: Easting (US survey foot)
- Y[north]: Northing (US survey foot)
Area of Use:
- name: United States (USA) - California - counties Kern; Los Angeles; San Bernardino; San Luis Obispo; Santa Barbara; Ventura.
- bounds: (-121.42, 32.76, -114.12, 35.81)
Coordinate Operation:
- name: SPCS83 California zone 5 (US survey foot)
- method: Lambert Conic Conformal (2SP)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
In [11]:
# Plot your parcel data with the column "LandValue"
ax = parcels.plot(column="LandValue", legend=True, figsize=FIGSIZE)

# Add a scale bar
plt.title("Land Value")
ax.add_artist(ScaleBar(1, scale_formatter = lambda value, unit: f"{value}000 ft"))
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[11], line 6
      4 # Add a scale bar
      5 plt.title("Land Value")
----> 6 ax.add_artist(ScaleBar(1, scale_formatter = lambda value, unit: f"{value}000 ft"))

NameError: name 'ScaleBar' is not defined
In [8]:
# Plot with a column "NetSecVal", Net Secured Value for each parcel.
parcels.plot(column="NetSecVal", legend = True, figsize=FIGSIZE)
Out[8]:
<Axes: >
In [9]:
# Visualize another column, "YearBuilt" and check out a problem
parcels.plot(column="YearBuilt", legend = True, figsize=FIGSIZE)
Out[9]:
<Axes: >
In [10]:
# Yes, there are two problems.
    # 1. Legend is too long. 
    # 2. Some parcels are disappeared.

# Therefore, let's see what's happening in the YearBuilt column.
# Firstly, let's print out the raw values.
parcels["YearBuilt"].values
Out[10]:
array(['1967', '2006', '1954', ..., '1973', '1973', '1961'], dtype=object)
In [11]:
# I see. They are in string type, not in numeric type. 
    # That's why the legend was categorical, not a continuous colorbar.

# So we are going to replace the YearBuilt columns, with the same values in float data type.
# You can convert the data type with the method, .astype()
parcels["YearBuilt"] = parcels["YearBuilt"].astype("float")
In [107]:
res_parcels['YearBuilt'].value_counts()[0]
Out[107]:
0
In [12]:
# Also, as there were some missing parcels, let's check if there's any Null values in YearBuilt column.
parcels["YearBuilt"].isna().sum()
Out[12]:
135
In [13]:
# Oh.. there were 135 null values, that's why there are some blanks. 
    # This is a data limitation or there must be some reason of the null values. 
    # So let's accept the fact that there exist null values.

# Now, plot the parcels GeoDataFrame with the "YearBuilt" column.
# Is there still something non-plausible..?
ax = parcels.plot(column="YearBuilt", legend = True, figsize=FIGSIZE)
plt.title("YearBuilt")
ax.add_artist(ScaleBar(1, scale_formatter = lambda value, unit: f"{value}000 ft"))
Out[13]:
<matplotlib_scalebar.scalebar.ScaleBar at 0x1d90250ed70>
In [14]:
# How can there be some buildings built so long ago..? Year 0..? 2024 years ago?!
# Let's check out the min and max values of the YearBuilt column values.
parcels["YearBuilt"].max(), parcels["YearBuilt"].min()
Out[14]:
(2015.0, 0.0)
In [15]:
# Certainly, it must mean that the parcels with the year value 0 are actually Null values. 
# So let's replace them with Null, so the plot will show better spatial variation.
parcels.loc[parcels["YearBuilt"]==0, "YearBuilt"] = None
In [101]:
### Your Turn!! Let's visualize the YearBuilt column with ScaleBar, once again. Does it look plausible?

ax = parcels.plot(column="YearBuilt", cmap = "viridis", legend = True, figsize=FIGSIZE)
plt.title("Year Built of Santa Barbara Land Parcels")
ax.add_artist(ScaleBar(1, scale_formatter = lambda value, unit: f"{value}000 ft"))
#plt.axis("off")
Out[101]:
<matplotlib_scalebar.scalebar.ScaleBar at 0x1d949f3a500>
In [124]:
import geopandas as gpd
import folium

# Load your GeoDataFrame and convert it to EPSG 4326
parcels = parcels.to_crs(4326)

# Check the occurrence of null values in the "YearBuilt" column
null_count = parcels['YearBuilt'].isnull().sum()
total_count = len(parcels)
print("Number of null values in 'YearBuilt':", null_count)
print("Total number of parcels:", total_count)

# Create GeoDataFrame for parcels with null and non-null values in "YearBuilt" column
parcels_null = parcels[parcels['YearBuilt'].isnull()]
parcels_not_null = parcels[~parcels['YearBuilt'].isnull()]

# Create a folium map
map_center = [parcels.unary_union.centroid.y, parcels.unary_union.centroid.x]
m = folium.Map(location=map_center, zoom_start=15, control_scale=True)

# Add GeoDataFrames to the map with different colors
folium.GeoJson(parcels_null, style_function=lambda x: {'fillColor': 'red', 'color': 'black'}).add_to(m)
folium.GeoJson(parcels_not_null, style_function=lambda x: {'fillColor': 'blue', 'color': 'black'}).add_to(m)

# Display the map
m
Number of null values in 'YearBuilt': 136
Total number of parcels: 1872
Out[124]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Step 2. Subset the Residential Parcels for Optimization¶

In [17]:
# Let's check what the dataset has in "LandUse" column using the .value_counts() method.
    # It returns the unique values and frequency.
parcels["LandUse"].value_counts()
Out[17]:
SINGLE FAMILY RESIDENCE              1560
CONDOS,COMMUNITY APT PROJS            182
RESIDENTIAL INCOME, 2-4 UNITS          42
APARTMENTS, 5 OR MORE UNITS            29
VACANT                                 25
PARKS                                   5
OFFICE BUILDINGS, MULTI-STORY           5
PARKING LOTS                            3
SUPERMARKETS                            2
RESTAURANTS,BARS                        1
SERVICE STATIONS                        1
RETAIL STORES, SINGLE STORY             1
CHURCHES, RECTORY                       1
STORE AND OFFICE COMBINATION            1
MIXED USE-COMMERCIAL/RESIDENTIAL        1
PROFESSIONAL BUILDINGS                  1
OFFICE BUILDINGS, SINGLE STORY          1
COMMERCIAL AND OFFICE CONDOS,PUDS       1
Name: LandUse, dtype: int64
In [18]:
# Let's take a subset of the original parcels data, sothat only residential parcels can be included.

# In .loc[], the first argument selects rows based on a condition, and the optional second argument specifies columns.
# The .isin() function generates a binary series indicating whether each value is present in the given array.
res_parcels = parcels.loc[parcels["LandUse"].isin(["SINGLE FAMILY RESIDENCE", "CONDOS,COMMUNITY APT PROJS", "RESIDENTIAL INCOME, 2-4 UNITS", "APARTMENTS, 5 OR MORE UNITS"])]
res_parcels
Out[18]:
OBJECTID APN LAYER Situs1 Situs2 Acreage LandUse UseCode PctTransf LandValue StrImpr LivImpr NetSecVal Net_Impr YearBuilt Bedrooms Bathrooms geometry
0 131101 041-350-024 Ground 2315 EDGEWATER WAY SANTA BARBARA, CA 93109 0.52 SINGLE FAMILY RESIDENCE 0100 100.0 288614 331840 0 620454 331840 1967.0 4 3.00 POLYGON ((6039997.884 1971525.657, 6039982.106...
1 131102 041-350-012 Ground 2307 EDGEWATER WAY SANTA BARBARA, CA 93109 0.18 SINGLE FAMILY RESIDENCE 0100 26.6 550902 1085141 0 1636043 1085141 2006.0 3 3.50 POLYGON ((6040044.305 1971477.086, 6040032.277...
2 131103 041-350-015 Ground 2211 EDGEWATER WAY SANTA BARBARA, CA 93109 0.59 SINGLE FAMILY RESIDENCE 0100 100.0 1641000 175000 0 1816000 175000 1954.0 2 2.50 POLYGON ((6040237.820 1971288.705, 6040179.849...
3 131104 041-350-016 Ground 2201 EDGEWATER WAY SANTA BARBARA, CA 93109 0.50 SINGLE FAMILY RESIDENCE 0100 100.0 1644611 1482404 0 3127015 1482404 1948.0 3 3.00 POLYGON ((6040320.008 1971271.353, 6040268.735...
4 131105 041-350-017 Ground 2111 EDGEWATER WAY SANTA BARBARA, CA 93109 0.76 SINGLE FAMILY RESIDENCE 0100 0.0 63587 156682 0 213269 156682 1950.0 4 3.00 POLYGON ((6040450.139 1971243.880, 6040398.978...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
1867 250105 045-290-008 Condo Third Floor 30 BARRANCA AVE 3 SANTA BARBARA, CA 93109 0.00 CONDOS,COMMUNITY APT PROJS 0323 100.0 46816 70592 0 117408 70592 1973.0 2 1.75 POLYGON ((6048328.471 1972754.300, 6048325.112...
1868 250106 045-290-009 Condo Third Floor 36 BARRANCA AVE 6 SANTA BARBARA, CA 93109 0.00 CONDOS,COMMUNITY APT PROJS 0323 100.0 750000 668000 0 1418000 668000 1973.0 2 1.75 POLYGON ((6048287.470 1972741.338, 6048291.580...
1869 250107 045-290-010 Condo Third Floor 36 BARRANCA AVE 5 SANTA BARBARA, CA 93109 0.00 CONDOS,COMMUNITY APT PROJS 0323 100.0 309553 206366 0 515919 206366 1973.0 2 2.00 POLYGON ((6048283.551 1972700.429, 6048271.984...
1870 250108 045-290-011 Condo Third Floor 40 BARRANCA AVE 3 SANTA BARBARA, CA 93109 0.00 CONDOS,COMMUNITY APT PROJS 0323 100.0 124884 187149 0 312033 187149 1973.0 2 1.75 POLYGON ((6048234.146 1972677.397, 6048238.386...
1871 250109 045-330-008 Condo Third Floor 101 OCEANO AVE 5 SANTA BARBARA, CA 93109 0.00 CONDOS,COMMUNITY APT PROJS 0332 100.0 305105 298615 0 603720 298615 1961.0 2 1.75 POLYGON ((6048280.057 1972978.876, 6048296.690...

1813 rows × 18 columns

In [21]:
# Also, let's say that we want the parcels where YearBuilt is identified.
# Check out how many parcels don't have YearBuilt identified.
res_parcels["YearBuilt"].isna().sum()
Out[21]:
0
In [22]:
# Then, let's drop them.
res_parcels = res_parcels.loc[-res_parcels["YearBuilt"].isna()]
res_parcels = res_parcels.reset_index(drop=True)
res_parcels
Out[22]:
OBJECTID APN LAYER Situs1 Situs2 Acreage LandUse UseCode PctTransf LandValue StrImpr LivImpr NetSecVal Net_Impr YearBuilt Bedrooms Bathrooms geometry
0 131101 041-350-024 Ground 2315 EDGEWATER WAY SANTA BARBARA, CA 93109 0.52 SINGLE FAMILY RESIDENCE 0100 100.0 288614 331840 0 620454 331840 1967.0 4 3.00 POLYGON ((6039997.884 1971525.657, 6039982.106...
1 131102 041-350-012 Ground 2307 EDGEWATER WAY SANTA BARBARA, CA 93109 0.18 SINGLE FAMILY RESIDENCE 0100 26.6 550902 1085141 0 1636043 1085141 2006.0 3 3.50 POLYGON ((6040044.305 1971477.086, 6040032.277...
2 131103 041-350-015 Ground 2211 EDGEWATER WAY SANTA BARBARA, CA 93109 0.59 SINGLE FAMILY RESIDENCE 0100 100.0 1641000 175000 0 1816000 175000 1954.0 2 2.50 POLYGON ((6040237.820 1971288.705, 6040179.849...
3 131104 041-350-016 Ground 2201 EDGEWATER WAY SANTA BARBARA, CA 93109 0.50 SINGLE FAMILY RESIDENCE 0100 100.0 1644611 1482404 0 3127015 1482404 1948.0 3 3.00 POLYGON ((6040320.008 1971271.353, 6040268.735...
4 131105 041-350-017 Ground 2111 EDGEWATER WAY SANTA BARBARA, CA 93109 0.76 SINGLE FAMILY RESIDENCE 0100 0.0 63587 156682 0 213269 156682 1950.0 4 3.00 POLYGON ((6040450.139 1971243.880, 6040398.978...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
1731 250105 045-290-008 Condo Third Floor 30 BARRANCA AVE 3 SANTA BARBARA, CA 93109 0.00 CONDOS,COMMUNITY APT PROJS 0323 100.0 46816 70592 0 117408 70592 1973.0 2 1.75 POLYGON ((6048328.471 1972754.300, 6048325.112...
1732 250106 045-290-009 Condo Third Floor 36 BARRANCA AVE 6 SANTA BARBARA, CA 93109 0.00 CONDOS,COMMUNITY APT PROJS 0323 100.0 750000 668000 0 1418000 668000 1973.0 2 1.75 POLYGON ((6048287.470 1972741.338, 6048291.580...
1733 250107 045-290-010 Condo Third Floor 36 BARRANCA AVE 5 SANTA BARBARA, CA 93109 0.00 CONDOS,COMMUNITY APT PROJS 0323 100.0 309553 206366 0 515919 206366 1973.0 2 2.00 POLYGON ((6048283.551 1972700.429, 6048271.984...
1734 250108 045-290-011 Condo Third Floor 40 BARRANCA AVE 3 SANTA BARBARA, CA 93109 0.00 CONDOS,COMMUNITY APT PROJS 0323 100.0 124884 187149 0 312033 187149 1973.0 2 1.75 POLYGON ((6048234.146 1972677.397, 6048238.386...
1735 250109 045-330-008 Condo Third Floor 101 OCEANO AVE 5 SANTA BARBARA, CA 93109 0.00 CONDOS,COMMUNITY APT PROJS 0332 100.0 305105 298615 0 603720 298615 1961.0 2 1.75 POLYGON ((6048280.057 1972978.876, 6048296.690...

1736 rows × 18 columns

In [106]:
res_parcels['YearBuilt'].isnull().sum()
Out[106]:
0
In [37]:
# Check the output GeoDataFrame, anything non-plausible?
res_parcels["Bedrooms"].value_counts()
Out[37]:
3     904
2     401
4     338
5      54
1      23
6      11
7       2
24      1
0       1
8       1
Name: Bedrooms, dtype: int64
In [ ]:
# How would you fix this? 
In [39]:
# Then, with the residential parcels having YearBuilt identified,
# let's check the spatial distribution of the number of Bedrooms. Plot with a column "Bedrooms"

res_parcels.plot("Bedrooms", figsize=(15,5), legend=True, categorical=True, cmap="YlOrRd")
plt.axis("off")
plt.title("The number of Bedrooms")
Out[39]:
Text(0.5, 1.0, 'The number of Bedrooms')
In [40]:
# Let's find the parcels, which doesn't have any Bedroom.
res_parcels.loc[res_parcels.Bedrooms==0]
Out[40]:
OBJECTID APN LAYER Situs1 Situs2 Acreage LandUse UseCode PctTransf LandValue StrImpr LivImpr NetSecVal Net_Impr YearBuilt Bedrooms Bathrooms geometry
1139 155252 045-185-006 Ground 1419 SHORELINE DR SANTA BARBARA, CA 93109 0.32 SINGLE FAMILY RESIDENCE 0100 100.0 77851 128819 0 199670 128819 1946.0 0 2.0 POLYGON ((6044967.574 1971082.045, 6044995.993...
In [43]:
# Let's create a new column, "bedroom" in the GeoDataFrame. (column name is case-sensitive)
    # The "bedroom" column will have a binary (True/False) value indicating if the parcel has greater than zero Bedrooms.
res_parcels["bedroom"] = res_parcels.Bedrooms > 0
    
# Now visualize the binary column, "bedroom".
res_parcels.plot("bedroom", figsize=(15,5), legend=True, categorical=True, cmap="YlOrRd")
plt.axis("off")
plt.title("parcels w/ > 1 bedroom")
Out[43]:
Text(0.5, 1.0, 'parcels w/ > 1 bedroom')
In [102]:
res_parcels.loc[res_parcels.Bedrooms==0]
Out[102]:
OBJECTID APN LAYER Situs1 Situs2 Acreage LandUse UseCode PctTransf LandValue StrImpr LivImpr NetSecVal Net_Impr YearBuilt Bedrooms Bathrooms geometry bedroom selected
1139 155252 045-185-006 Ground 1419 SHORELINE DR SANTA BARBARA, CA 93109 0.32 SINGLE FAMILY RESIDENCE 0100 100.0 77851 128819 0 199670 128819 1946.0 0 2.0 POLYGON ((-119.71290 34.39651, -119.71280 34.3... False False

Step 3. Optimization using GurobiPy: Knapsack Problem¶

Knapsack Problem Equations¶

Notations:

The objective is to maximize the total value of items selected while keeping the total weight within the knapsack capacity.

Let:

  • ( $n$ ) be the number of items.
  • ( $v_i$ ) be the value of the ( i )th item.
  • ( $w_i$ ) be the weight of the ( i )th item.
  • ( $W$ ) be the maximum capacity of the knapsack.

where:

  • ( $x_i$ ) is a binary decision variable:
    • ( $x_i$ = 1 ) if the ( i )th item is selected.
    • ( $x_i$ = 0 ) otherwise.

Complete Formulation:

The complete formulation of the Knapsack problem is as follows:

Maximize:

$ \sum_{i=1}^{n} v_i \cdot x_i $

Subject to:

$ \sum_{i=1}^{n} w_i \cdot x_i \leq W $

$ x_i \in \{0, 1\}, \quad \forall i = 1, 2, \ldots, n $

In [29]:
pip install gurobipy
Requirement already satisfied: gurobipy in /Users/ezrabanks/anaconda3/lib/python3.11/site-packages (10.0.3)
Note: you may need to restart the kernel to use updated packages.
In [30]:
# Let's import the GurobiPy into this Notebook.

import gurobipy as gp
from gurobipy import GRB
In [105]:
# Define the data (weights, values, and knapsack capacity)
objective_weights = res_parcels['Bedrooms']  # Example weights of items
constraint_weights = res_parcels['LandValue']   # Example values of items
BUDGET = 1000000            # Knapsack capacity

# Create a new Gurobi model
m = gp.Model("knapsack")

# Add decision variables for each item (binary variables indicating whether to include the item in the knapsack or not)
x = m.addVars(len(res_parcels), vtype=GRB.BINARY, name="x")

# Set objective function: maximize total value
m.setObjective(sum(objective_weights[i] * x[i] for i in range(len(res_parcels))), sense=GRB.MAXIMIZE)

# Add constraint: total weight of selected items must not exceed knapsack capacity (BUDGET)
#m.addConstr(sum(constraint_weights[i] * x[i] for i in range(len(res_parcels))) <= BUDGET)

# Optimize the model
m.optimize()
Gurobi Optimizer version 11.0.1 build v11.0.1rc0 (win64 - Windows 10.0 (19045.2))

CPU model: 12th Gen Intel(R) Core(TM) i7-12700, instruction set [SSE2|AVX|AVX2]
Thread count: 12 physical cores, 20 logical processors, using up to 20 threads

Optimize a model with 0 rows, 1736 columns and 0 nonzeros
Model fingerprint: 0x39170ef2
Variable types: 0 continuous, 1736 integer (1736 binary)
Coefficient statistics:
  Matrix range     [0e+00, 0e+00]
  Objective range  [1e+00, 2e+01]
  Bounds range     [1e+00, 1e+00]
  RHS range        [0e+00, 0e+00]
Found heuristic solution: objective 5271.0000000

Explored 0 nodes (0 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 1 (of 20 available processors)

Solution count 1: 5271 

Optimal solution found (tolerance 1.00e-04)
Best objective 5.271000000000e+03, best bound 5.271000000000e+03, gap 0.0000%
In [48]:
# Now, add the model result to the GeoDataFrame
res_parcels["selected"] = [x[i].x > 0.5 for i in range(len(res_parcels))]

# Then plot the result
res_parcels.plot("selected", figsize=(15,5), legend=True, cmap="rainbow", edgecolor="black")
plt.axis("off")
plt.title("Selected Residential Parcels")
Out[48]:
Text(0.5, 1.0, 'Selected Residential Parcels')
In [50]:
! pip install folium
Defaulting to user installation because normal site-packages is not writeable
Collecting folium
  Downloading folium-0.16.0-py2.py3-none-any.whl (100 kB)
     -------------------------------------- 100.0/100.0 kB 2.9 MB/s eta 0:00:00
Requirement already satisfied: requests in c:\tools\anaconda3\lib\site-packages (from folium) (2.28.1)
Collecting xyzservices
  Downloading xyzservices-2024.4.0-py3-none-any.whl (81 kB)
     ---------------------------------------- 82.0/82.0 kB ? eta 0:00:00
Collecting branca>=0.6.0
  Downloading branca-0.7.1-py3-none-any.whl (25 kB)
Requirement already satisfied: numpy in c:\tools\anaconda3\lib\site-packages (from folium) (1.23.5)
Requirement already satisfied: jinja2>=2.9 in c:\tools\anaconda3\lib\site-packages (from folium) (3.1.2)
Requirement already satisfied: MarkupSafe>=2.0 in c:\tools\anaconda3\lib\site-packages (from jinja2>=2.9->folium) (2.1.1)
Requirement already satisfied: certifi>=2017.4.17 in c:\tools\anaconda3\lib\site-packages (from requests->folium) (2024.2.2)
Requirement already satisfied: charset-normalizer<3,>=2 in c:\tools\anaconda3\lib\site-packages (from requests->folium) (2.0.4)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in c:\tools\anaconda3\lib\site-packages (from requests->folium) (1.26.14)
Requirement already satisfied: idna<4,>=2.5 in c:\tools\anaconda3\lib\site-packages (from requests->folium) (3.4)
Installing collected packages: xyzservices, branca, folium
Successfully installed branca-0.7.1 folium-0.16.0 xyzservices-2024.4.0
In [111]:
import folium

res_parcels = res_parcels.to_crs(4326)
# Create a Folium map centered around the mean latitude and longitude of the GeoDataFrame
map_center = [res_parcels.unary_union.centroid.y, res_parcels.unary_union.centroid.x]
m = folium.Map(location=map_center, zoom_start=15, control_scale=True)


# Define a style function to set colors based on the binary column
def style_function(feature):
    color = 'red' if feature['properties']['selected'] == 1 else 'blue'
    return {'fillColor': color, 'color': color, 'weight': 1, 'fillOpacity': 0.6}

tooltip = folium.GeoJsonTooltip(
    fields=["APN", "selected", "LandValue", "Bedrooms"]
)


# Add GeoDataFrame to the map with the defined style and tooltip functions
folium.GeoJson(
    res_parcels,
    style_function=style_function,
    tooltip=tooltip
).add_to(m)

# Display the map
m
Out[111]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [ ]:
res_parcels = res_parcels.to_crs(4326)
# Create a Folium map centered around the mean latitude and longitude of the GeoDataFrame
map_center = [res_parcels.unary_union.centroid.y, res_parcels.unary_union.centroid.x]
m = folium.Map(location=map_center, zoom_start=15, control_scale=True)

def style_function(feature):
    color = 'red' if feature['properties']['selected'] == 1 else 'blue'
    return {'fillColor': color, 'color': color, 'weight': 1, 'fillOpacity': 0.6}

Chapter 2. Wildfires and Major Roads¶

Step 1. Download the Road Network Data in Santa Barbara¶

Note: OSMnx libarary provides the filtering functionality when downloading a graph.
However, this time we will not use it to practice data processing in Geopandas.

Hint: Santa Barbara Bounding Coordinates

- outhernmost Latitude: 34.2714° N <br>
- Northernmost Latitude: 35.5295° N <br>
- Westernmost Longitude: -120.6976° W <br>
- Easternmost Longitude: -118.9600° W <br>
In [53]:
# Fetch the road network data into a variable, "road"
north, south, east, west = 35.5295, 34.2714, -118.9600, -120.6976
graph = ox.graph_from_bbox(north, south, east, west, 
                           network_type='drive')
#                           custom_filter='["highway"~"primary|trunk|motorway"]')

# Plot the road network
ox.plot_graph(graph)

# Convert it into GeoDataFrame
C:\Users\Student\AppData\Local\Temp\ipykernel_8252\862118477.py:3: FutureWarning: The `north`, `south`, `east`, and `west` parameters are deprecated and will be removed in the v2.0.0 release. Use the `bbox` parameter instead. See the OSMnx v2 migration guide: https://github.com/gboeing/osmnx/issues/1123
  graph = ox.graph_from_bbox(north, south, east, west,
Out[53]:
(<Figure size 800x800 with 1 Axes>, <Axes: >)
In [55]:
# Convert it into GeoDataFrame
roads = ox.graph_to_gdfs(graph, nodes=False, edges=True)
In [56]:
roads.plot()
Out[56]:
<Axes: >
In [57]:
# Check the unique values in "highway" column to select only Major Roads.
roads['highway'].value_counts()
Out[57]:
residential                               90457
tertiary                                  13278
secondary                                  7567
unclassified                               4555
primary                                    3027
motorway_link                               813
trunk                                       626
motorway                                    590
secondary_link                              371
primary_link                                156
[residential, unclassified]                 146
living_street                               136
[residential, tertiary]                     113
tertiary_link                               112
trunk_link                                   92
[motorway, trunk]                            17
[tertiary, unclassified]                     14
[secondary, tertiary]                        11
[residential, living_street]                  8
road                                          8
[secondary, unclassified]                     4
[motorway, primary]                           4
crossing                                      4
[secondary, secondary_link]                   3
escape                                        3
[residential, secondary, unclassified]        2
disused                                       2
[residential, secondary]                      2
[residential, motorway_link]                  2
[tertiary_link, tertiary]                     2
[tertiary_link, residential]                  2
[residential, tertiary, unclassified]         2
[motorway_link, primary]                      1
[motorway_link, secondary_link]               1
[motorway, motorway_link, primary]            1
[trunk_link, primary_link]                    1
[motorway_link, trunk_link]                   1
Name: highway, dtype: int64
In [58]:
# Try to plot the road GeoDataFrame with a column "highway"
roads.plot["highway"]
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[58], line 2
      1 # Try to plot the road GeoDataFrame with a column "highway"
----> 2 roads.plot["highway"]

TypeError: 'GeoplotAccessor' object is not subscriptable
In [59]:
# Error? No worries. That's intended.
# Why do you think error happens? Hint: read the error message.
roads["highway_"] = roads["highway"].astype("str")
In [60]:
# How can we solve this issue? Let's stringify the column, "highway" and create a new column, "highway_"


# Then plot with the new column, "highway_"
roads.plot("highway_", legend=True, figsize=(20,20))
plt.axis("off")
Out[60]:
(-120.78512326500001, -118.872378035, 34.20852533, 35.59235327)
In [63]:
# Now, let's create another column, "major", indicating whether each road segment is part of Major Roads or not.
    # Let's consider the "primary" or "trunk" or "motorway" as Major Roads
roads["major"] = roads["highway_"].str.contains("primary|trunk|motorway")
    
# Then plot with the new column, "major"
roads.plot("major",  legend=True, figsize=(20,20), cmap="Spectral_r")
plt.axis("off")
Out[63]:
(-120.78512326500001, -118.872378035, 34.20852533, 35.59235327)
In [65]:
# Create a subset GeoDataFrame, "major_roads" by selecting the rows where major column value is True

major_roads = roads.loc[roads["major"]==True]
major_roads
major_roads.crs
Out[65]:
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
In [66]:
# Interactive plot with the major roads. Do you agree that those are major roads?
# Interactive plot with the major roads. Do you agree that those are major roads?

import folium

# Create a Folium map centered around the mean latitude and longitude of the GeoDataFrame
map_center = [major_roads.unary_union.centroid.y, major_roads.unary_union.centroid.x]
m = folium.Map(location=map_center, zoom_start=7, control_scale=True)


# Define a style function to set colors based on the binary column
def style_function(feature):
    color = 'red' if feature['properties']['selected'] == 1 else 'blue'
    return {'fillColor': color, 'color': color, 'weight': 1, 'fillOpacity': 0.6}

tooltip = folium.GeoJsonTooltip(
    fields=["highway_"]
)


# Add GeoDataFrame to the map with the defined style and tooltip functions
folium.GeoJson(
    major_roads,
    tooltip=tooltip,
    name="Major Roads"
).add_to(m)

folium.LayerControl().add_to(m)
# Display the map
m
Out[66]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Step 2. Relate Ignition locations with Major Roads¶

In [82]:
# Load a shapefile, "Lab1Data/SBC_2007_2021_Wildland_Ignitions.shp" into a variable, "ignitions"
ignitions = gpd.read_file("Lab1Data/SBC_2007_2021_Wildland_Ignitions.shp")
In [81]:
# Plot the major_roads and ignitions together
In [ ]:
# Why it doesn't work? What should you check?


# Yes. Coordinate System. 
In [83]:
# Match the major_roads crs with the ignitions crs.
major_roads = major_roads.to_crs(ignitions.crs)
In [84]:
# Create a new column, "mr_distance" in ignitions GeoDataFrame. 
    # The column contains the distance from each ignition point to the closest major_roads segment.
    # Divide the distance values with 5280 (feet -> miles)
ignitions["mr_distance"] = ignitions.distance(major_roads.unary_union)
ignitions["mr_distance"] = ignitions["mr_distance"]/5280 # feet to miles
A = ignitions.plot("mr_distance", legend=True)
major_roads.plot(ax=A)
    
# Plot the result!
Out[84]:
<Axes: >
In [85]:
# Now, let's draw a histogram showing the distribution of Distance.
plt.figure(figsize=(10, 6))
ignitions.mr_distance.hist(bins=40, color='skyblue', edgecolor='black', alpha=0.7)

# Add titles and labels
plt.title('Histogram of Ignition to Major Road Distance')
plt.xlabel('Distance to Major Road (miles)')
plt.ylabel('Frequency')

# Add grid
plt.grid(axis='y', linestyle='--', alpha=0.7)

# Add a background color
plt.gca().set_facecolor('#f5f5f5')

# Add a legend (optional)
plt.legend(['Distance from Ignition to Major Road'])

# Show plot
plt.show()
In [86]:
# It's hard to see the result because of an outlier in Channel Island.
# What will be a good way to drop it?

ignitions = ignitions.drop(ignitions.mr_distance.argmax())
In [87]:
# After dropping the outlier, create a histogram once again.

plt.figure(figsize=(10, 6))
ignitions.mr_distance.hist(bins=40, color='skyblue', edgecolor='black', alpha=0.7)

plt.title('Histogram of Ignition to Major Road Distance')
plt.xlabel('Distance to Major Road (miles)')
plt.ylabel('Frequency')

plt.grid(axis='y', linestyle='--', alpha=0.7)

plt.gca().set_facecolor('#f5f5f5')

plt.legend(['Distance from Ignition to Major Road'])

plt.show()
In [88]:
# This time, visualize the ignition points with the distance value, and also plot major_roads together.
A = ignitions.plot("mr_distance",legend=True, figsize=(10,10), cmap="RdYlGn_r", legend_kwds={'shrink': 0.3})
major_roads.plot(ax=A, color="darkgrey", lw=0.5)
plt.axis("off")
plt.title("Distance from Ignition to Major Road (miles)")
Out[88]:
Text(0.5, 1.0, 'Distance from Ignition to Major Road (miles)')

Step 3. Area Burnt and Distance to Major Roads¶

Also, let's check out if the distance to Major Road is related with the area burnt by the fire.
What would be a good way to do so?

In [89]:
# Check out the columns contained in ignitions GeoDataFrame.
# Which column is having the area burnt information?
ignitions
Out[89]:
Id Alarm_Date Alarm_Time Cause Contain_Da Contain_Ti Incident_L Incident_N Acres POINT_X POINT_Y DPA FED_Cause Incident_1 geometry mr_distance
0 1 20070102 1513 7 20070102 1538 2825 Baseline Ave CA-SBC-0000044 2.00 -120.103645 34.638207 1 7 None POINT (5928955.913 2061259.921) 0.877190
1 2 20070412 1709 23 20070412 1732 3849 Roblar Ave CA-SBC-0002850 0.50 -120.069939 34.653624 1 2 None POINT (5939207.727 2066659.037) 1.305005
2 3 20070411 0843 5 20070411 0955 1135 Jonata Park Rd CA-SBC-0002803 0.50 -120.192122 34.633199 1 5 None POINT (5902309.367 2060006.431) 0.263935
3 4 20070123 1622 99 20070123 1651 6888 W Hwy 246 CA-SBC-0000689 0.10 -120.290366 34.633472 1 9 None POINT (5872766.228 2060764.421) 0.049065
4 5 20070427 0413 23 20070427 0617 E Hwy 246 @ Drum Canyon Rd CA-SBC-0003267 0.50 -120.288363 34.632188 1 2 None POINT (5873357.831 2060283.745) 0.034620
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
1208 0 10/27/2021 5:09:00 PM 17:09 99 10/27/2021 17:31 Hwy 1 CA-SBC-13553 0.10 -120.277002 34.522425 1 9 Cruces POINT (5875873.534 2020268.950) 0.007770
1209 0 11/4/2021 2:16:00 PM 14:16 2 11/4/2021 15:00 1470 N San Marcos Rd CA-SBC-13901 0.10 -119.801320 34.466034 1 2 San Marcos POINT (6018756.283 1996843.809) 1.367655
1210 0 11/14/2021 9:41:00 AM 9:41 99 11/14/2021 10:10 4100 Harris Grade CA-SBC-14367 0.17 -120.440616 34.735279 1 9 Harris POINT (5828485.635 2098866.021) 2.338022
1211 0 11/21/2021 6:28:00 PM 18:28 99 11/21/2021 21:00 5225 Figueroa Mtn Rd CA-SBC-14727 0.75 -120.087634 34.745155 1 9 Ranch POINT (5934582.955 2100070.912) 5.308591
1212 0 5/20/2021 9:02:00 PM 21:02 7 5/23/2021 None Loma Alta Dr CA-STB-5950 5.60 -119.706602 34.411130 2 7 Loma POINT (6046962.114 1976367.883) 0.223830

1212 rows × 16 columns

In [90]:
# Draw scatterplot with the mr_distance and area burnt columns.
plt.scatter(ignitions["mr_distance"], ignitions["Acres"])
Out[90]:
<matplotlib.collections.PathCollection at 0x1d9b32f3820>
In [91]:
# Using a threshold, select a subset of ignitions to exclude outliers.

ignitions_wo_outliers = ignitions.loc[ignitions["Acres"] < 50000]
In [92]:
# Draw scatterplot with the mr_distance and area burnt columns.
plt.figure(figsize=FIGSIZE)
plt.scatter(ignitions_wo_outliers["mr_distance"], ignitions_wo_outliers["Acres"])
plt.xlabel("Distance to Major Road")
plt.ylabel("Area Burnt in Acres")
plt.title("Scatter plot")
Out[92]:
Text(0.5, 1.0, 'Scatter plot')
In [93]:
# This time, let's use another way to drop the outliers. Use the IQR based approach.
# You can ask chatgpt for the function.

import numpy as np

def remove_outliers_iqr(data):
    """
    Remove outliers from a dataset using the Interquartile Range (IQR) method.
    
    Parameters:
    - data: A list or NumPy array containing the dataset.
    
    Returns:
    - Cleaned dataset with outliers removed.
    """
    # Convert data to NumPy array for easier manipulation
    data = np.array(data)
    
    # Calculate the first quartile (Q1)
    Q1 = np.percentile(data, 25)
    
    # Calculate the third quartile (Q3)
    Q3 = np.percentile(data, 75)
    
    # Calculate the interquartile range (IQR)
    IQR = Q3 - Q1
    
    # Define the lower bound and upper bound for outlier detection
    lower_bound = Q1 - 1.5 * IQR
    upper_bound = Q3 + 1.5 * IQR
    
    # Filter out the outliers
    cleaned_index = (data >= lower_bound) & (data <= upper_bound)
    
    return cleaned_index


ignitions_wo_outliers = ignitions.loc[remove_outliers_iqr(ignitions["Acres"])]
In [94]:
# Visualzie the scatter plot again. Don't forget to add titles and labels.

plt.figure(figsize=FIGSIZE)
plt.scatter(ignitions_wo_outliers["mr_distance"], ignitions_wo_outliers["Acres"])
plt.xlabel("Distance to Major Road")
plt.ylabel("Area Burnt in Acres")
Out[94]:
Text(0, 0.5, 'Area Burnt in Acres')
In [27]:
iv_parcels = gpd.read_file("Lab1Data/Parcel_IV.shp")
iv_parcels
Out[27]:
OBJECTID APN LAYER Situs1 Situs2 Acreage LandUse UseCode PctTransf LandValue StrImpr LivImpr NetSecVal Net_Impr YearBuilt Bedrooms Bathrooms geometry
0 130716 075-010-022 Ground 6875 EL COLEGIO RD GOLETA, CA 93117 9.60 SCHOOLS 7200 0.0 0 0 0 0 0 None 0 0.0 POLYGON ((5997851.219 1979476.854, 5998251.155...
1 130717 075-010-039 Ground None None 0.83 VACANT 0090 100.0 0 0 0 0 0 None 0 0.0 POLYGON ((5998411.116 1979466.012, 5998407.771...
2 130718 075-010-040 Ground None None 0.35 VACANT 0090 100.0 0 0 0 0 0 None 0 0.0 POLYGON ((5998631.052 1979460.936, 5998407.771...
3 130719 075-010-037 Ground None None 21.60 VACANT 0090 100.0 0 0 0 0 0 None 0 0.0 POLYGON ((5998749.878 1979458.194, 5998727.559...
4 130720 075-010-028 Ground 6795 EL COLEGIO RD GOLETA, CA 93117 11.78 VACANT 0090 0.0 0 0 0 0 0 None 0 0.0 POLYGON ((5999395.829 1979443.282, 5999388.511...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
912 244502 075-230-001 Condo First Floor 6636 DEL PLAYA DR # 1 GOLETA, CA 93117 0.03 CONDOS,COMMUNITY APT PROJS 0300 100.0 308718 308718 0 617436 308718 2003 1 1.0 POLYGON ((6000659.415 1976779.651, 6000658.790...
913 244503 075-230-003 Condo First Floor 6636 DEL PLAYA DR # 3 GOLETA, CA 93117 0.03 CONDOS,COMMUNITY APT PROJS 0300 100.0 308718 308718 0 617436 308718 2003 1 1.0 POLYGON ((6000659.596 1976790.249, 6000659.432...
914 248854 075-230-002 Condo Second Floor 6636 DEL PLAYA DR # 2 GOLETA, CA 93117 0.03 CONDOS,COMMUNITY APT PROJS 0300 100.0 308718 308718 0 617436 308718 2003 1 2.0 POLYGON ((6000646.610 1976774.894, 6000646.553...
915 248855 075-230-004 Condo Second Floor 6636 DEL PLAYA DR # 4 GOLETA, CA 93117 0.03 CONDOS,COMMUNITY APT PROJS 0300 100.0 308718 308718 0 617436 308718 2003 1 2.0 POLYGON ((6000659.898 1976810.348, 6000659.596...
916 239415 075-020-039 Right of Way (Private) None None 0.52 HIGHWAYS AND STREETS 8800 100.0 154 0 0 154 0 None 0 0.0 POLYGON ((6000356.910 1978984.122, 6000353.435...

917 rows × 18 columns

In [28]:
ax = iv_parcels.plot(column="NetSecVal", legend=True, figsize=FIGSIZE)

plt.title("NetSecVal")
ax.add_artist(ScaleBar(1, scale_formatter = lambda value, unit: f"{value}000 ft"))
Out[28]:
<matplotlib_scalebar.scalebar.ScaleBar at 0x13570c0d0>

Final Report ( Total 30 pts )¶

Chapter 1.¶

  1. Submit your final map visualization for the "YearBuilt" column.
    Ensure your map includes necessary elements such as title, legend, and scale bar.
    Do not turn off the axis.
    Utilize a suitable colormap. ( 3 pts )

  2. Load a new GeoDataFrame, "Lab1Data/Parcel_IV.shp".
    Visualize its NetSecValue column with Scalebar.
    You shouldn't turn off the axis for this map. ( 4 pts )

  3. Report the number of SB parcels selected in the given Knapsack Problem, aimed at maximizing the number of bedrooms within a budget of $1,000,000. ( 1 pt )

  4. Report the number of bedrooms that can be obtained in the given Knapsack Problem, with the budget of $1,000,000. ( 1 pt )

  5. Discuss the presence of null values in the "YearBuilt" column.
    Hint: Utilize the Folium library to overlay parcels with null or not-null values on a basemap. ( 2 pts )

  6. In a knapsack problem, what happens if the objective function sets "MINIMIZE" rather than "MAXIMIZE"? ( 2 pt )

  7. How many bathrooms does the parcel with zero bedrooms have? ( 1 pt )

  8. In a knapsack problem, what happens if there's no constraint? Hint: Comment out the setConstraint row ( 2 pts )

  9. Evaluate whether the current knapsack approach is appropriate or not, and provide reasoning for your conclusion. ( 2 pts )

  10. Design a new Knapsack problem scenario, referring to the description.txt file for details on parcel data columns.
    m. vbnm Explain the scenario you have created. ( 2 pts )

  11. Run the Knapsack Problem code with your custom scenario.
    Visualize the results in an interactive plot using Folium.
    Ensure that the color scheme is different from blue/red.
    (Optional) You can try the different parcels in Isla Vista region (parcl_IV.shp) instead of the parcels used for instruction. ( 4 pts )

Chapter 2.¶

  1. Submit a screenshot of your histogram and scatter plot without outliers.
    Your plots must have title, x and y labels. (4 pts)
  1. What insight can you get out of the map, histogram and scatterplot you created in Chapter 2?
    How can you utilize the insights for the urban wildfire risk management? ( 2 pts )
In [53]:
# Define the data (weights, values, and knapsack capacity)
objective_weights = iv_parcels['Acreage']  # Example weights of items
constraint_weights = iv_parcels['NetSecVal']   # Example values of items
BUDGET = 1000000            # Knapsack capacity

# Create a new Gurobi model
m = gp.Model("knapsack")

# Add decision variables for each item (binary variables indicating whether to include the item in the knapsack or not)
x = m.addVars(len(iv_parcels), vtype=GRB.BINARY, name="x")

# Set objective function: maximize total value
m.setObjective(sum(objective_weights[i] * x[i] for i in range(len(iv_parcels))), sense=GRB.MAXIMIZE)

# Add constraint: total weight of selected items must not exceed knapsack capacity (BUDGET)
m.addConstr(sum(constraint_weights[i] * x[i] for i in range(len(iv_parcels))) <= BUDGET)

# Optimize the model
m.optimize()
Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (mac64[x86])

CPU model: Intel(R) Core(TM) i5-5350U CPU @ 1.80GHz
Thread count: 2 physical cores, 4 logical processors, using up to 4 threads

Optimize a model with 1 rows, 917 columns and 840 nonzeros
Model fingerprint: 0x18b6a785
Variable types: 0 continuous, 917 integer (917 binary)
Coefficient statistics:
  Matrix range     [2e+02, 1e+08]
  Objective range  [3e-02, 2e+01]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+06, 1e+06]
Found heuristic solution: objective 86.3530000
Presolve removed 0 rows and 531 columns
Presolve time: 0.01s
Presolved: 1 rows, 386 columns, 386 nonzeros
Variable types: 0 continuous, 386 integer (371 binary)
Found heuristic solution: objective 88.7830000

Root relaxation: objective 9.088949e+01, 1 iterations, 0.00 seconds (0.00 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0   90.88949    0    1   88.78300   90.88949  2.37%     -    0s
H    0     0                      90.0930000   90.88949  0.88%     -    0s
H    0     0                      90.7230000   90.88949  0.18%     -    0s
H    0     0                      90.8330000   90.88949  0.06%     -    0s
     0     0   90.87853    0    1   90.83300   90.87853  0.05%     -    0s

Cutting planes:
  Cover: 1

Explored 1 nodes (2 simplex iterations) in 0.17 seconds (0.01 work units)
Thread count was 4 (of 4 available processors)

Solution count 5: 90.833 90.723 90.093 ... 86.353

Optimal solution found (tolerance 1.00e-04)
Best objective 9.083300000000e+01, best bound 9.083300000000e+01, gap 0.0000%
In [54]:
# Now, add the model result to the GeoDataFrame
iv_parcels["selected"] = [x[i].x > 0.5 for i in range(len(iv_parcels))]

# Then plot the result
iv_parcels.plot("selected", figsize=(15,5), legend=True, cmap="plasma", edgecolor="black")
plt.axis("off")
plt.title("Selected Residential Parcels")
Out[54]:
Text(0.5, 1.0, 'Selected Residential Parcels')
In [55]:
import folium

iv_parcels = iv_parcels.to_crs(4326)
# Create a Folium map centered around the mean latitude and longitude of the GeoDataFrame
map_center = [iv_parcels.unary_union.centroid.y, iv_parcels.unary_union.centroid.x]
m = folium.Map(location=map_center, zoom_start=15, control_scale=True)


# Define a style function to set colors based on the binary column
def style_function(feature):
    color = 'yellow' if feature['properties']['selected'] == 1 else 'blue'
    return {'fillColor': color, 'color': color, 'weight': 1, 'fillOpacity': 0.6}

tooltip = folium.GeoJsonTooltip(
    fields=["APN", "selected", "LandValue", "Acreage"]
)


# Add GeoDataFrame to the map with the defined style and tooltip functions
folium.GeoJson(
    iv_parcels,
    style_function=style_function,
    tooltip=tooltip
).add_to(m)

# Display the map
m
Out[55]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [ ]:
 
In [ ]: